<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>Ensembles of oscillators</title>
<link rel="stylesheet" href="../../../../../../../doc/src/boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../../index.html" title="Chapter 1. Boost.Numeric.Odeint">
<link rel="up" href="../tutorial.html" title="Tutorial">
<link rel="prev" href="lattice_systems.html" title="Lattice systems">
<link rel="next" href="using_boost__units.html" title="Using boost::units">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../logo.jpg"></td>
<td align="center"><a href="../../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="lattice_systems.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="using_boost__units.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h3 class="title">
<a name="boost_numeric_odeint.tutorial.ensembles_of_oscillators"></a><a class="link" href="ensembles_of_oscillators.html" title="Ensembles of oscillators">Ensembles
      of oscillators</a>
</h3></div></div></div>
<p>
        Another important high dimensional system of coupled ordinary differential
        equations is an ensemble of <span class="emphasis"><em>N</em></span> all-to-all coupled phase
        oscillators <a class="link" href="../literature.html#synchronization_pikovsky_rosenblum">[9] </a>.
        It is defined as
      </p>
<p>
        <span class="emphasis"><em>dφ<sub>​k</sub> / dt = ω<sub>​k</sub> + ε / N Σ<sub>​j</sub> sin( φ<sub>​j</sub> - φ<sub>​k</sub> )</em></span>
      </p>
<p>
        The natural frequencies <span class="emphasis"><em>ω<sub>​i</sub></em></span> of each oscillator follow
        some distribution and <span class="emphasis"><em>ε</em></span> is the coupling strength. We
        choose here a Lorentzian distribution for <span class="emphasis"><em>ω<sub>​i</sub></em></span>. Interestingly
        a phase transition can be observed if the coupling strength exceeds a critical
        value. Above this value synchronization sets in and some of the oscillators
        oscillate with the same frequency despite their different natural frequencies.
        The transition is also called Kuramoto transition. Its behavior can be analyzed
        by employing the mean field of the phase
      </p>
<p>
        <span class="emphasis"><em>Z = K e<sup>i Θ</sup> = 1 / N Σ<sub>​k</sub>e<sup>i φ<sub>​k</sub></sup></em></span>
      </p>
<p>
        The definition of the system function is now a bit more complex since we
        also need to store the individual frequencies of each oscillator.
      </p>
<p>
</p>
<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">vector</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">&gt;</span> <span class="identifier">container_type</span><span class="special">;</span>


<span class="identifier">pair</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">&gt;</span> <span class="identifier">calc_mean_field</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">container_type</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">)</span>
<span class="special">{</span>
    <span class="identifier">size_t</span> <span class="identifier">n</span> <span class="special">=</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">size</span><span class="special">();</span>
    <span class="keyword">double</span> <span class="identifier">cos_sum</span> <span class="special">=</span> <span class="number">0.0</span> <span class="special">,</span> <span class="identifier">sin_sum</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span>
    <span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">i</span><span class="special">&lt;</span><span class="identifier">n</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span>
    <span class="special">{</span>
        <span class="identifier">cos_sum</span> <span class="special">+=</span> <span class="identifier">cos</span><span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">);</span>
        <span class="identifier">sin_sum</span> <span class="special">+=</span> <span class="identifier">sin</span><span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">);</span>
    <span class="special">}</span>
    <span class="identifier">cos_sum</span> <span class="special">/=</span> <span class="keyword">double</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">);</span>
    <span class="identifier">sin_sum</span> <span class="special">/=</span> <span class="keyword">double</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">);</span>

    <span class="keyword">double</span> <span class="identifier">K</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span> <span class="identifier">cos_sum</span> <span class="special">*</span> <span class="identifier">cos_sum</span> <span class="special">+</span> <span class="identifier">sin_sum</span> <span class="special">*</span> <span class="identifier">sin_sum</span> <span class="special">);</span>
    <span class="keyword">double</span> <span class="identifier">Theta</span> <span class="special">=</span> <span class="identifier">atan2</span><span class="special">(</span> <span class="identifier">sin_sum</span> <span class="special">,</span> <span class="identifier">cos_sum</span> <span class="special">);</span>

    <span class="keyword">return</span> <span class="identifier">make_pair</span><span class="special">(</span> <span class="identifier">K</span> <span class="special">,</span> <span class="identifier">Theta</span> <span class="special">);</span>
<span class="special">}</span>


<span class="keyword">struct</span> <span class="identifier">phase_ensemble</span>
<span class="special">{</span>
    <span class="identifier">container_type</span> <span class="identifier">m_omega</span><span class="special">;</span>
    <span class="keyword">double</span> <span class="identifier">m_epsilon</span><span class="special">;</span>

    <span class="identifier">phase_ensemble</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">n</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">g</span> <span class="special">=</span> <span class="number">1.0</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">epsilon</span> <span class="special">=</span> <span class="number">1.0</span> <span class="special">)</span>
    <span class="special">:</span> <span class="identifier">m_omega</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_epsilon</span><span class="special">(</span> <span class="identifier">epsilon</span> <span class="special">)</span>
    <span class="special">{</span>
        <span class="identifier">create_frequencies</span><span class="special">(</span> <span class="identifier">g</span> <span class="special">);</span>
    <span class="special">}</span>

    <span class="keyword">void</span> <span class="identifier">create_frequencies</span><span class="special">(</span> <span class="keyword">double</span> <span class="identifier">g</span> <span class="special">)</span>
    <span class="special">{</span>
        <span class="identifier">boost</span><span class="special">::</span><span class="identifier">mt19937</span> <span class="identifier">rng</span><span class="special">;</span>
        <span class="identifier">boost</span><span class="special">::</span><span class="identifier">cauchy_distribution</span><span class="special">&lt;&gt;</span> <span class="identifier">cauchy</span><span class="special">(</span> <span class="number">0.0</span> <span class="special">,</span> <span class="identifier">g</span> <span class="special">);</span>
        <span class="identifier">boost</span><span class="special">::</span><span class="identifier">variate_generator</span><span class="special">&lt;</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">mt19937</span><span class="special">&amp;,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">cauchy_distribution</span><span class="special">&lt;&gt;</span> <span class="special">&gt;</span> <span class="identifier">gen</span><span class="special">(</span> <span class="identifier">rng</span> <span class="special">,</span> <span class="identifier">cauchy</span> <span class="special">);</span>
        <span class="identifier">generate</span><span class="special">(</span> <span class="identifier">m_omega</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">m_omega</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">gen</span> <span class="special">);</span>
    <span class="special">}</span>

    <span class="keyword">void</span> <span class="identifier">set_epsilon</span><span class="special">(</span> <span class="keyword">double</span> <span class="identifier">epsilon</span> <span class="special">)</span> <span class="special">{</span> <span class="identifier">m_epsilon</span> <span class="special">=</span> <span class="identifier">epsilon</span><span class="special">;</span> <span class="special">}</span>

    <span class="keyword">double</span> <span class="identifier">get_epsilon</span><span class="special">(</span> <span class="keyword">void</span> <span class="special">)</span> <span class="keyword">const</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">m_epsilon</span><span class="special">;</span> <span class="special">}</span>

    <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">container_type</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">container_type</span> <span class="special">&amp;</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="comment">/* t */</span> <span class="special">)</span> <span class="keyword">const</span>
    <span class="special">{</span>
        <span class="identifier">pair</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">&gt;</span> <span class="identifier">mean</span> <span class="special">=</span> <span class="identifier">calc_mean_field</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">);</span>
        <span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">i</span><span class="special">&lt;</span><span class="identifier">x</span><span class="special">.</span><span class="identifier">size</span><span class="special">()</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span>
            <span class="identifier">dxdt</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">m_omega</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">m_epsilon</span> <span class="special">*</span> <span class="identifier">mean</span><span class="special">.</span><span class="identifier">first</span> <span class="special">*</span> <span class="identifier">sin</span><span class="special">(</span> <span class="identifier">mean</span><span class="special">.</span><span class="identifier">second</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">);</span>
    <span class="special">}</span>
<span class="special">};</span>
</pre>
<p>
      </p>
<p>
        Note, that we have used <span class="emphasis"><em>Z</em></span> to simplify the equations
        of motion. Next, we create an observer which computes the value of <span class="emphasis"><em>Z</em></span>
        and we record <span class="emphasis"><em>Z</em></span> for different values of <span class="emphasis"><em>ε</em></span>.
      </p>
<p>
</p>
<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">statistics_observer</span>
<span class="special">{</span>
    <span class="keyword">double</span> <span class="identifier">m_K_mean</span><span class="special">;</span>
    <span class="identifier">size_t</span> <span class="identifier">m_count</span><span class="special">;</span>

    <span class="identifier">statistics_observer</span><span class="special">(</span> <span class="keyword">void</span> <span class="special">)</span>
    <span class="special">:</span> <span class="identifier">m_K_mean</span><span class="special">(</span> <span class="number">0.0</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_count</span><span class="special">(</span> <span class="number">0</span> <span class="special">)</span> <span class="special">{</span> <span class="special">}</span>

    <span class="keyword">template</span><span class="special">&lt;</span> <span class="keyword">class</span> <span class="identifier">State</span> <span class="special">&gt;</span>
    <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">State</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span>
    <span class="special">{</span>
        <span class="identifier">pair</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">&gt;</span> <span class="identifier">mean</span> <span class="special">=</span> <span class="identifier">calc_mean_field</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">);</span>
        <span class="identifier">m_K_mean</span> <span class="special">+=</span> <span class="identifier">mean</span><span class="special">.</span><span class="identifier">first</span><span class="special">;</span>
        <span class="special">++</span><span class="identifier">m_count</span><span class="special">;</span>
    <span class="special">}</span>

    <span class="keyword">double</span> <span class="identifier">get_K_mean</span><span class="special">(</span> <span class="keyword">void</span> <span class="special">)</span> <span class="keyword">const</span> <span class="special">{</span> <span class="keyword">return</span> <span class="special">(</span> <span class="identifier">m_count</span> <span class="special">!=</span> <span class="number">0</span> <span class="special">)</span> <span class="special">?</span> <span class="identifier">m_K_mean</span> <span class="special">/</span> <span class="keyword">double</span><span class="special">(</span> <span class="identifier">m_count</span> <span class="special">)</span> <span class="special">:</span> <span class="number">0.0</span> <span class="special">;</span> <span class="special">}</span>

    <span class="keyword">void</span> <span class="identifier">reset</span><span class="special">(</span> <span class="keyword">void</span> <span class="special">)</span> <span class="special">{</span> <span class="identifier">m_K_mean</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span> <span class="identifier">m_count</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> <span class="special">}</span>
<span class="special">};</span>
</pre>
<p>
      </p>
<p>
        Now, we do several integrations for different values of <span class="emphasis"><em>ε</em></span>
        and record <span class="emphasis"><em>Z</em></span>. The result nicely confirms the analytical
        result of the phase transition, i.e. in our example the standard deviation
        of the Lorentzian is 1 such that the transition will be observed at <span class="emphasis"><em>ε =
        2</em></span>.
      </p>
<p>
</p>
<pre class="programlisting"><span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">n</span> <span class="special">=</span> <span class="number">16384</span><span class="special">;</span>
<span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">dt</span> <span class="special">=</span> <span class="number">0.1</span><span class="special">;</span>

<span class="identifier">container_type</span> <span class="identifier">x</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">);</span>

<span class="identifier">boost</span><span class="special">::</span><span class="identifier">mt19937</span> <span class="identifier">rng</span><span class="special">;</span>
<span class="identifier">boost</span><span class="special">::</span><span class="identifier">uniform_real</span><span class="special">&lt;&gt;</span> <span class="identifier">unif</span><span class="special">(</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">2.0</span> <span class="special">*</span> <span class="identifier">M_PI</span> <span class="special">);</span>
<span class="identifier">boost</span><span class="special">::</span><span class="identifier">variate_generator</span><span class="special">&lt;</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">mt19937</span><span class="special">&amp;,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uniform_real</span><span class="special">&lt;&gt;</span> <span class="special">&gt;</span> <span class="identifier">gen</span><span class="special">(</span> <span class="identifier">rng</span> <span class="special">,</span> <span class="identifier">unif</span> <span class="special">);</span>

<span class="comment">// gamma = 1, the phase transition occurs at epsilon = 2</span>
<span class="identifier">phase_ensemble</span> <span class="identifier">ensemble</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">,</span> <span class="number">1.0</span> <span class="special">);</span>
<span class="identifier">statistics_observer</span> <span class="identifier">obs</span><span class="special">;</span>

<span class="keyword">for</span><span class="special">(</span> <span class="keyword">double</span> <span class="identifier">epsilon</span> <span class="special">=</span> <span class="number">0.0</span> <span class="special">;</span> <span class="identifier">epsilon</span> <span class="special">&lt;</span> <span class="number">5.0</span> <span class="special">;</span> <span class="identifier">epsilon</span> <span class="special">+=</span> <span class="number">0.1</span> <span class="special">)</span>
<span class="special">{</span>
    <span class="identifier">ensemble</span><span class="special">.</span><span class="identifier">set_epsilon</span><span class="special">(</span> <span class="identifier">epsilon</span> <span class="special">);</span>
    <span class="identifier">obs</span><span class="special">.</span><span class="identifier">reset</span><span class="special">();</span>

    <span class="comment">// start with random initial conditions</span>
    <span class="identifier">generate</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">gen</span> <span class="special">);</span>

    <span class="comment">// calculate some transients steps</span>
    <span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">runge_kutta4</span><span class="special">&lt;</span> <span class="identifier">container_type</span> <span class="special">&gt;()</span> <span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">ref</span><span class="special">(</span> <span class="identifier">ensemble</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">10.0</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">);</span>

    <span class="comment">// integrate and compute the statistics</span>
    <span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">runge_kutta4</span><span class="special">&lt;</span> <span class="identifier">container_type</span> <span class="special">&gt;()</span> <span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">ref</span><span class="special">(</span> <span class="identifier">ensemble</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">100.0</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">ref</span><span class="special">(</span> <span class="identifier">obs</span> <span class="special">)</span> <span class="special">);</span>
    <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">epsilon</span> <span class="special">&lt;&lt;</span> <span class="string">"\t"</span> <span class="special">&lt;&lt;</span> <span class="identifier">obs</span><span class="special">.</span><span class="identifier">get_K_mean</span><span class="special">()</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="special">}</span>
</pre>
<p>
      </p>
<p>
        The full cpp file for this example can be found here <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/phase_oscillator_ensemble.cpp" target="_top">phase_oscillator_ensemble.cpp</a>
      </p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2009-2015 Karsten Ahnert and Mario Mulansky<p>
        Distributed under the Boost Software License, Version 1.0. (See accompanying
        file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
      </p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="lattice_systems.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="using_boost__units.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
